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NOMENCLATURE 


a,b,k,h 

A 


A,B,C1,D,E,F,H 

C 

CURV 

e 

GRAD 

H.O.T. 

i 

o 

n 

JW 

K 

K 

£ 

L 

n 

N 

NN 

'’a 

Q 

r 

K 

e 

RMS 

t 

At 

II 

X 

Ax 

©,@.0.  etc. 


Variables  in  two-dimensional  Taylor  series  expansion 
method 

Amplification  factor  in  linear  stability  method,  A 
=  modulus;  coefficient  in  QUICK  method  of  Leonard 
(1979) 

Coefficients  in  implicit,  finite-difference  methods 

Courant  number 

Curvature 

A  subscript  meaning  exact  equation 
Gradient 

Higher  order  truncation  error  terms 
Space  step  index 
Imaginary  number  index 

One  half  the  base  width  of  the  test  triangle 

Wave  number,  27t/L 

Diffusion  (dispersion)  coefficient 

Left  side  wall  value  as  subscript 

Wave  length 

Time  step  index  as  superscript;  also  number  of  Taylor 
series  terms 

Total  number  of  grid  points  in  one  full  wave  length 
Total  number  of  time  steps 
Peclet  number 
Phase  ratio 

Right  side  wall  value  as  subscript 
Exact  continuum  response  function 
Root-mean-square  value 

Time;  as  subscript  means  time  differentiation 

Time  step 

Velocity 

Space  coordinate;  as  subscript  means  space 
differentiation 

Space  step 

Truncation  error  terms 
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One - d imens i ona 1 
Two-dimensional 
Dimensionless  wave  number 

Dimensionless  diffusion  coefficient  (the  diffusion 
number) 

Weighting  factor,  0  <  6  <  1  in  implicit  schemes 
Dependent  variables 
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AN  IMPLICIT,  WIGGLE-FREE,  AND  ACCURATE  UPSTREAM 
FINITE-DIFFERENCE  ALGORITHM  FOR  THE 


ONE-DIMENSIONAL  TRANSPORT-DIFFUSION  EQUATION 

PART  I ;  INTRODUCTION 

1.  Implicit  finite-difference  algorithms  use  more  than  one  un¬ 
known  variable  at  the  next  time  level  (n+1)  in  the  calculation.  They 
thereby  require  the  simultaneous  solution  of  many  equations  at  time 
(n+1)  in  order  to  advance  the  computation  for  marching-type,  initial 
value/boundary  value  problems. 

2.  In  general,  implicit  schemes  have  better  stability  properties 
than  explicit  schemes.  For  propagation-dominated  problems,  this  means 
running  with  Courant  numbers  greater  than  unity  to  use  more  economic 
time  steps  in  the  computation.  In  general,  the  price  paid  is  less  accu¬ 
racy  and  larger  Central  Processing  Unit  (CPU)  times  on  the  computer.  In 
addition,  implicit  schemes  require  a  prescribed  order  or  direction  of 
the  computation  in  the  space  of  the  independent  variables.  This  can  be 
called  the  algorithmic  structure  of  the  scheme  (Abbott  1979).  Use  of  an 
inappropriate  structure  will  result  in  an  ill-posed  problem  and  unstable 
solution . 

3.  Eng"  ers  require  robust  schemes  in  their  computational  sys¬ 
tems  to  pr  greater  flexibility  in  use,  yet  maintain  acceptable 

standards  oi  ..cxuracy  over  a  wide  variety  of  conditions.  Properly  de¬ 
signed  implicit  schemes  have  met  these  needs  in  many  aspects  of  Compu¬ 
tational  Hydraulics. 

4.  For  transport-diffusion  computations,  Leonard  (1979)  has  de¬ 
veloped  the  explicit"  QUICKEST  modeling  procedure.  It  avoids  the  **wiggle” 
instability  problem  associated  with  central  differencing  of  the  advec- 
tion  term.  It  also  eliminates  the  inaccuracies  of  numerical  diffusion 
resulting  when  only  first-order,  upstream  differencing  procedures  are  em¬ 
ployed.  The  purpose  of  this  research  effort  was  to  develop  and  study 
numerical  properties  of  an  implicit  algorithm  consistent  with  the  Leonard 
(1979)  scheme  for  the  one-dimensional  transport-diffusion  equation. 
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5.  Part  II  of  this  report  reviews  the  explicit  scheme  of  Leonard 
(1979)  and  proves  its  equivalence  to  first-order  differences  with 
judicious  removal  of  the  truncation  error  terms.  In  Part  III,  two  im¬ 
plicit  schemes  are  developed  and  their  characteristics  are  investigated 
using  linear  stability  analysis  methods  and  standard  numerical  tests. 
Both  implicit  methods  are  then  compared  with  the  explicit  Leonard 
scheme.  Conclusions  and  recommendations  follow  in  Part  IV. 
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PART  II:  EXPLICIT  LEONARD  ALGORITHM  (QUICKEST) 


6.  It  is  highly  instructive  to  briefly  review  Leonardos  (1979) 
explicit  scheme  for  unsteady  transport  diffusion,  which  is  labeled 
QUICKEST  (quadratic  upstream  interpolation  for  convective  kinematics 
with  estimated  streaming  terms).  Consider  the  differential  equation 

^t  ^  ^  ^^xx  constants  (1) 

together  with  appropriate  initial  and  boundary  conditions.*  A  scalar 
quantity  is  converted  by  constant  velocity  u  and  diffused  by  a 
constant  diffusion  coefficient  K  in  one  space  dimension.  A  source 
term  could  also  be  added,  if  necessary. 

7.  For  the  finite-difference  grid  shown  in  Figure  1,  Leonard 

(1979)  proposed  using  a  basic,  three-point,  upstream-weighted,  quadratic 
interpolation  scheme  for  "wall'*  values  (right)  and  (t)^^  (left)  of 

a  control  volume  when  u  is  positive  to  the  right.  For  steady-state 
conditions,  Leonard  took 


Figure  1.  Schematic  of  QUICKEST  algorithm  for  flow  from 

left  to  right 


*  See  "Nomenclature"  on  page  5  for  definitions  of  symbols  and 
abbreviations . 
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♦r  =  ♦i»l/2  =  *  ♦.)  -  -  24,.  t  ♦..j)  (2! 

♦i  =  ^i-i/a  =  *  ♦i-i>  -  *<♦!  -  ^♦i.,  *  ♦i-j)  (3) 

with  A  =  1/8  for  steady  flows  (QUICK  algorithm) . 

8.  For  unsteady-state  conditions,  Equations  2  and  3  were  incor¬ 
porated  in  an  exact  integral  formulation  where  average  wall  values  over 
time  increment  At  were  employed.  Complete  details  are  beyond  the 
intended  scope  of  this  review. 


Explicit  Operator  (QUICKEST) 


9.  The  resultant  explicit,  finite-difference  operator  for  the 
QUICKEST  algorithm  becomes  (Leonard  1979,  p.  80,  Equation  57) 


1  1 


r[i/2(«“^l  +  ^  C-GRAD^  -  ^  (1  -  -  anCURV^j  ^ 

\^[l/2(fl  +  CGRAD^  -  ^  (1  .  C2  -  SDCURV^jj  (4) 

l^^AxGRAD^  -  ^  C-CURVj.^  -  ^AxGRAD^j  -  ~  C'CURVj^J 


where ; 


C  =  -7 —  =  the  Courant  Number 


r  =  — 2  ”  dimensionless  diffusion  coefficient,  or  the 
Ax^  diffusion  number 


GRAD  = 
r 


^1+1 


CURV  = 
r 


,  -  24?  +  4?., 
•^i+l  ^i+1 


for  flow  from  left  to  right  for  both  walls  (r  =  right;  £  =  left).  Note 
that  for  this  case,  both  curvature  terms  CURV  are  actually  centered 
upstream  by  one-half  increment  and  not  at  r  or  £  as  indicated.  This 
is  a  key  aspect  of  the  QUICKEST  scheme.  In  Equation  4,  the  physical 
diffusion  is  calculated  using  a  standard,  centered,  second  difference 
operator  at  point  i  .  The  overall  (global)  truncation  error  is  third 
order  in  space,  and  the  scheme  gives  a  conservative  formulation  with  no 
wiggle  instabilities  present. 


Equivalence  to  a  Higher  Order  Accuracy  Method 


10.  As  shall  be  demonstrated  below,  QUICKEST  is  equivalent  to  a 
forward-time,  centered-space  (FTCS)  scheme  in  which  the  truncation 
errors  resulting  from  all  terms  are  subtracted  out  in  a  prescribed 
manner . 

11.  From  a  Taylor series  expansion  about  point  (i,n)  in  Fig¬ 
ure  1,  one  obtains  the  finite-difference  analog  to  Equation  1. 
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I 


,1 


@ 


-  K 


^xxxx 


+  ♦  + 

6 !  ^xxxxxx 


H.O.T.] 


(5) 


A  forward-time  difference  (term  (^)  and  a  centered-space  difference 
(term  are  employed  along  with  a  centered,  second  difference  (term 

for  the  diffusion.  The  bracketed  terms  and  are 

the  truncation  errors  associated  with  terms 
respectively. 

12.  To  eliminate  the  time  derivatives  in  Equation  1  can  be 

further  differentiated  in  time  and  space  to  yield 


(h  =  -  2uK^ 

^tt  ^xx  ^xxx 


A  =  -u^6 
^ttt  ^xxx 


(6) 

(7) 


and 


■^tx 


^'xtt 


u(|»  -f  K6 
^xx  ^xxx 

(a) 

=  U  0 

^xxx 

(b) 

(8) 


when  4th  derivatives  are  neglected.  Equation  8  will  be  used  later  in 
this  report.  Substituting  Equations  6  and  7  into  5  and  grouping  like- 
ordered  differentials  gives 


- at—  ^  - 2^ - 


Cl  - 


i+1 


Ax' 


*  #  [♦xj" 


4! 


A  + 
^xxxx 


-  ^)[*xx.]” 
H.O.T.] 


(9) 


The  last  bracketed  term  in  Equation  9  is  the  truncation  error 
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3  2 

0(At  ,  Ax  )  and  contains  fourth-order  and  higher  differentials  which 
are  neglected.  Multiplying  through  by  At  and  clearing  all  n-level 
terms  to  the  right-hand  side  (RHS)  gives 


n 

i-1 


+  C 


'XX 


(10) 


The  bracketed  term  in  Equation  10  is  the  remaining  truncation  error 
below  3nd  the  method  in  which  it  is  finite-differenced  is  cen¬ 

tral  to  the  equivalency  of  Equations  4  and  10. 

13.  Using  the  terms  defined  by  Leonard  (1979)  in  Equation  4 

gives 


) 

^xx 


GRAD  -  GRAD, 
r  Z 

Ax 


2(t»"  + 
1 


Ax 


n 

i-1 


or 


^(''*xx)i  = 

Thus,  for  ,  the  centering  is  precisely  at  (i,n)  as  required  by 

Equation  10,  Now,  if  the  centering  of  is  judiciously  moved  up¬ 

stream  to  (i  -  1/2,  n),  it  is  assumed  that 


II  I] 

V^xxx/  .  --  \  ^xxx/  . 


i-1/2 


CURV  -  CURV. 
r _ Z 

Ax 
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.  3^“  +  3^“ 


(W)"  (♦xxx)”^^^^  =  -  <^t 

The  identity 


(12) 


>'2  [♦"  ♦  Cl]-  [♦;  ^ 


(13) 


is  also  required. 

14.  Putting  Equations  11,  12,  and  13  into  Equation  10  gives 


.n+l 

"•i 


-  C  [l/2  (<!>;  +  -  1/2  ((|.“  +  J  ]  +  r[(  AxGRAdJ 

-  ^AxGRADj^^j+  C  1^1  Ax  (gRAD^  -  GRAD^^) 


+  ^  (l  -  -  6r 


)  (^CURV^  CURVj^  )]  (14) 


The  final  trick  to  recover  the  form  given  by  Leonard  is  to  split  the 
last  term  in  Equation  14  simply  by  taking 

-6r  =  -3r  -  3r 


and  through  some  further  rearrangement.  Equation  4  is  thus  recovered 
and  the  derivation  is  complete. 

15.  For  programming  purposes,  a  much  simpler  version  can  be  ob¬ 
tained  by  inserting  the  difference  Equations  11  and  12  for  and 

*^xxx  *  respectively.  After  grouping  all  like  terms,  one  obtains 


14 


(15) 


(■ 


2r  + 


+ 


c 

2 


(1 


§  (i  -  -  6r)j  4,"  ,  ,  r  *  ^ 

-  -  6r)]  (i  -  -  6r)]  4,“.2 


16.  The  algorithm  for  QUICKEST  does  reduce  to  that  given  by 
Leonard  (1979)  for  the  steady-state  algorithm  (QUICK)  when  term  for 
time  derivative  truncation  error  is  omitted.  Thus  QUICKEST  implicitedly 
takes  A  =  1/6  in  Equations  2  and  3.  Other  reduced  forms  are 

Pure  Advection  (F  =  0) 

=  I*"  -  I  (C^  -  3C  +  2)  +  I  (C^  -  2C  -  1)(|»“ 

-  I  (C^  -  C  -  2)  +  §  (C^  -  l)<l'i.2  (16) 

Pure  Diffusion  (C  =  0) 

♦r' '  <  *  r  (17) 

When  r  =  0  and  C  =  +1  ,  Equation  16  reduces  to 
point-to-point  transport  for  an  exact  result. 


Linear  Stability  Analysis 


17.  As  a  further  review  and  check  of  these  results,  a  von  Neumann 
linear  stability  analysis  was  made  of  Equation  15.  This  resulted  in  the 
following  expression  for  the  amplification  factor  A 


A(a)  =  1  +  2(r  +  i  C^)(cos  a  -  1) 


+  j^g  C(1  -  -  6r)(4  cos  Of  -  3  -  cos  2a)J 

-n  C  [sin  a  +  g  (1  -  -  6r)(2  sin  a  -  sin  2a )J 


(18) 
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where : 


a  =  kAx  *  the  dimensionless  wave  number 

Ott 

K  =  the  wave  number  — 

“  N 

N  =  the  number  of  grid  intervals  per  wave  length  L 

o 

n  =  the  imaginary  number  index 

Equation  18  is  identical  to  that  derived  by  Leonard  (1979,  p.  81, 
Equation  59),  which  gives  further  confidence  in  these  results. 

18.  An  example  amplitude  portrait  for  the  pure  advection  case  is 
shown  in  Figure  2  for  three  values  of  the  Courant  number.  The  modulus 
of  A  versus  a  plot  (Figure  2a)  gives  more  resolution  for  the  higher 
wave  numbers  (a  ■*  n)  and  consequently  has  been  employed  for  all  subse¬ 
quent  portraits.  No  phase  errors  exist  for  C  =  0.5  and  1.0  (F  =  0)  . 

19.  The  complete  stability  range  in  the  (r,C)  plane  for  the 

QUICKEST  algorithm  as  mapped  by  Leonard  (1979)  is  reproduced  in  Fig¬ 
ure  3.  Lines  of  constant  Peclet  niunber  defined  as 


(19) 


are  also  shown.  When  F  =  0  (P^  =  »)  ,  the  QUICKEST  scheme  is  unstable 
when  1  <  C  <  2  and  also  for  C  >  2  .  The  scheme  is  stable  for  all 
other  F,  C  values  within  the  space  shown  with  the  P^  lines.  Al¬ 
though  this  explicit  scheme  is  stable  for  a  considerable  range  above 
C  =  1  ,  it  also  is  unstable  in  a  "corner"  beyond  F  =0.5  and  below 
C  =  0.5  .  It  would  be  anticipated  that  an  implicit  version  would  elim¬ 
inate  this  unstable  corner. 

20.  But  Figure  3  reveals  nothing  about  the  accuracy  of  the  scheme. 
For  the  continuum  Equation  1,  the  exact  amplitude  response  function 
R^(a)  is  given  by 

Rg(o)  =  exp(-Fo(^)  (20) 

since  the  modulus  for  pure  advection  alone  is  unity  over  all  a  .  The 
exact  phase  ratio  Q^  is  also  unity  for  all  a  since  pure  diffusion  has 
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oo 


(b)  When  wave  number  space  Is  defined  by  N  . 

Figure  2.  Amplitude  portraits  of  QUICKEST  for  C  =  1.0, 
0.9,  and  0.5  for  pure  advection  (F  =  0) 
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no  phase  properties.  Examples  of  the  accuracy  characteristics  of  QUICK¬ 
EST  are  presented  in  Figure  4  when  C=0.5,r=l.l  and  in  Figure  5 
for  C=0.5  ,r=1.2.  The  exact  response  is  both  cases  is  also  pre¬ 
sented.  For  the  stable  case  (1.1,  0.5),  large  differences  occur  for 
both  1  A  I  and  Q  at  higher  wave  numbers  meaning  low  accuracy  for  these 
short  wave  components.  The  unstable  case  (1.2,  0.5)  is  clearly  re¬ 
vealed  by  the  amplitude  portrait  in  Figure  5. 

21.  It  would  be  very  beneficial  to  characterize  the  accuracy  of 
QUICKEST  by  one  number  over  the  entire  (r,C)  plane.  To  do  this,  the 
root-mean-square  (RMS)  value  of  the  differences  between  continuum 
(exact)  and  discrete  (approximate)  responses  was  computed  using  59 
values  equally  spaced  in  N  .  Iso-RMS  "difference”  lines  were  then  con¬ 
structed  to  give  the  "contour"  type  plot  for  amplitude  response  pre¬ 
sented  in  Figure  6.  The  weighted-average  RMS  difference  value  for  the 
stable  area  (also  traced  on  Figure  6)  is  0.0568.  In  N-space,  this  RMS 
norm  was  found  dependent  on  the  number  of  differences  employed.  Further 
research  is  needed  to  refine  this  approach. 

Verification 


22.  A  convenient  and  informative  numerical  test  for  verification 
of  transport-diffusion  schemes  is  the  propagation  of  test  triangles  with 
periodic  boundary  conditions.  Since  numerical  diffusion  and  convective 
phase  errors  are  of  primary  concern,  these  tests  are  made  with  no  physi¬ 
cal  diffusion  (F  =  0)  present.  Two  examples  are  shown  in  Figure  7. 
Triangles  with  peak  of  1.0  unit  centered  at  i  =  50  on  a  100-unit  grid 
and  half-base  widths  of  20  units  (Figure  7a)  and  2  units  (Figure  7b)  are 
propagated  with  C  =  +0.5  (left  to  right).  With  this  Courant  number 
and  using  periodic  boundary  conditions,  200  time  steps  are  needed  to 
return  the  triangle  to  its  initial  position  shown. 

23.  The  points  shown  in  Figures  7a  and  7b  are  results  using  the 
QUICKEST  algorithm  after  200  time  steps.  No  phase  errors  are  indicated, 
which  confirms  the  von  Neumann  analysis.  However,  numerical  amplitude 
errors  for  high  wave  number  components  cause  the  peak  to  be  dasiped. 
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Figure  4.  Amplitude  and  phase  portraits  of  QUICKEST  and 
QUICKIST  (implicit)  schemes  when  C=0.5,  r=l.l 


0 


2i 


Figure  5.  Amplitude  and  phase  portraits  of  QUICKEST  and 
QUICKIST  (implicit)  schemes  when  C=0.5  ,  r=  1.2  . 
Explicit  QUICKEST  is  unstable 
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Figure  6.  Contour  lines  of  equal  RMS  difference 
values  (iso-amplitudes)  on  the  (F,  C)  plane  for 
QUICKEST.  Stability  range  also  shown 


(a)  Wide  triangle  with  half-base  width  of  20  units. 


(b)  Narrow  triangle  with  half-base  width  of  2  units. 

Figure  7.  Propagation  test  results  for  pure 
advection  case  using  periodic  boundary  condi¬ 
tions  and  QUICKEST  algorithm  with  C  =  0.5 

especially  for  the  "spiked"  triangle  (JW  =  2)  example.  The  damping  of 
the  wider  triangle  (JW  =  20)  is  relatively  minor  when  compared  with  that 
produced  by  many  first-order  finite-difference  schemes  (see,  e.g.,  Reid 
and  fiasco  (1981)  for  other  examples). 

Wiggles 

24.  Wiggles  as  defined  by  Leonard  (1979)  are 

"...spatially  decaying  or  growing  oscillations  of 
wavelength  2Ax(i.e.  Of  =  7t,  N  =  2)  -  typical  of 
central  difference  solutions  of  the. . .convective- 
diffusion  problem  for  >  2  . 
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...Interestingly  enough,  the  criterion  for  the 
absence  of  wiggles,  ^  <  2  ,  is  not  a  stability 
condition  in  the  von  Newmann  sense....**  (p.  63) 

Roache  (1977)  discusses  this  as  the  **zero-overshoot**  requirement  while 
Abbott  (1979)  mentions  **zig-2agging**  of  free-surface  flow  schemes.  A 
full  discussion  of  the  sensitivity  of  the  advection  term  for  various 
schemes  and  sensitivity  requirements  to  prevent  wiggles  is  found  in 
Leonard  (1979).  Upstream  differencing  of  the  advection  term  always 
prevents  wiggles. 

25.  To  further  examine  the  wiggle  problem  for  centered  convection 

schemes,  an  implicit  FTCS  scheme  was  developed  with  no  attempt  to  remove 

the  truncation  errors.  An  implicit  scheme  is  necessary  for  stability 

since  the  explicit  FTCS  scheme  for  pure  advection  is  always  unstable. 

The  advection  term  is  centered  in  space  and  weighted  between  n  and 

n  +  1  time  levels  by  a  weighting  factor  0(0  <  0  <  1)  .  The  usual 

2  2 

case  employed  0  =  1/2  to  give  the  a(At  ,  Ax  )  for  the  truncation 
error. 

26.  A  linear  stability  analysis  proved  the  scheme  stable  for  all 
(r,  C)  values  when  0  >  0.5  .  Some  phase  errors  existed,  but  |  A  |  h  1 
for  all  or  when  C  =  0.5  and  1.0  (0  =  1/2)  .  Propagation  test  re¬ 
sults  for  the  wide  shape  are  shown  in  Figure  8.  Clearly,  an  oscilla¬ 
tion  of  length  2Ax  is  present  and  can  only  be  explained  as  a  **wiggle** 
instability  of  central  differencing  since  =  »  (F  =  0)  .  In  Fig¬ 
ure  8b,  the  Courant  number  is  1.0  requiring  only  100  time  steps  to  com¬ 
plete  the  propagation  cycle.  In  each  case,  the  wave  moves  too  slow  as 
predicted  by  the  phase  portrait  for  higher  wave  numbers.  Additional 
tests  with  ?^=  2  eliminated  the  wiggles.  When  P^  =  5  ,  wiggles 
were  barely  perceptible  near  the  peak  as  shown  in  Figure  9,  where  the 
triangles  only  represent  initial  conditions  since  physical  diffusion  is 
now  included. 
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0  10  20^^30  40  50  60  70  'SO'  'OO  ^^O 


(a)  When  C=0.5  and  after  2(X)  time  steps. 


0  10  20V30  40  50  60  70  80  90  100 


(b)  When  C=1.0  and  after  100  time  steps. 

Figure  8.  Example  of  wiggle  instability  for 
implicit,  centered  advection  of  test  triangle 
with  periodic  boundary  conditions 


(a)  When  C=0.5  and  after  200  time  steps. 


(b)  When  C=1.0  and  after  100  time  steps. 

Figure  9.  Suppression  of  wiggle  instability  when 
<  2  for  implicit,  centered  advection  of  test 
triangle.  When  *  5  ,  a  slight  oscillation  is 
barely  visible 


26 


PART  III:  IMPLICIT  ALGORITHMS 


27.  Many  alternative  ways  exist  to  formulate  implicit,  higher 
order  accurate  difference  schemes  for  Equation  1.  Two  methods  are  de¬ 
scribed  below,  both  of  which  utilize  elements  of  the  explicit,  Leonard 
(1979)  algorithm. 


Implicit  Leonard  Algorithm  for  Transport -Diffusion  (QUICKIST) 


28.  In  discussions  of  computational  efficiency  and  expense  of  the 

QUICK  (steady-state)  method,  Leonard  (1979,  p.  73)  states 

Implicit,  tridiagonal  procedures  are  similarly 
economical.  In  that  case,  one  merely  uses  (29),* 
treating  the  linear-interpolation  term  implicitly, 
with  the  curvature  evaluated  as  as  explicit 
source  term. 

As  derived  in  the  previous  chapter,  the  linear-interpolation  term  stems 
from  central  differencing  the  advection  term.  Consequently,  this  term 
is  made  implicit  through  use  of  a  weighting  coefficient  0  such  that 
0  <  0  <  1  ,  i.e. , 


n  ^i+1  " 

(1  -  =  (1  -  6) - J3;; - 


-  (1  -  ♦)  ^ * 


xxxxx  5 ! 


H.O.T. 

-*  i 


(21) 


.n+1  ^ 

.n+l  _  »  ***1+1  o 

®  ® 


[A  2  .4  T  n+1 

4'  +  4'  %■  +  H.O.T.  (22) 

^xxx  3!  ^xxxxx  5!  J  ^ 


A  similar  weighting  makes  the  diffusion  implicit  (Crank-Nicholson 
method)  so  that,  by  analogy  to  Equation  5, 


*  Equations  2  and  3  in  this  report. 
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©  ®  ,, 

r  r  2  ^n+1 

-K(l  -  9)  ^  V,.  ^  "  -  ''4i^  ♦.««]  i 

When  0  =  0,  terms  Q,  (?),  (?),  and  @  disappear  and  the  explicit 
scheme  given  by  Equation  5  is  recovered.  Again,  neglecting  all  fourth 
and  higher  order  derivatives.  Equation  23  gives,  after  rearrangement, 
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29.  The  numbered  brackets  identify  truncation  error  terms.  Trun¬ 

cation  error  time  derivatives  are  again  converted  to  space  derivatives 
using  Equations  6  and  7  to  give  terms  and  ^2^ ,  respectively.  All 

truncation  error  derivatives  are  at  time  level  n  except  term  A 

third  spatial  derivative  requires  four  grid  points ,  which  is  incompatible 
with  a  tridiagonal  algorithm.  Based  on  the  derivation  in  the  previous 
chapter  and  the  quote  above  (Leonard  1979,  p.  73),  for  evaluation  of  the 
curvature  as  an  "explicit  source  term,"  term  is  taken  at  time  level 

n  .  This  is  a  key  assumption  in  this  first  method.  Also,  the  time 
derivative  is  still  centered  at  grid  point  (i,n)  so  that  this  implicit 
method  reduces  to  the  explicit  method  of  Leonard  (1979)  called  QUICKEST 
when  0=0.  Here  the  underlined  E  also  could  stand  for  the  explicit 
algorithm.  By  analogy,  for  later  reference,  this  implicit  algorithm  is 
labeled  QUICKIST. 

30.  Removing  term  in  Equation  24  to  time  level  n  ,  multiply¬ 
ing  through  by  At  ,  and  rearranging  by  grouping  all  like  terms  gives 


y 
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Note  that  taking  term  to  level  n  effectively  removes  the  influence 
of  0  on  the  truncation  errors  associated  with  the  convection  term. 

31.  Finally,  if  again  the  assumption  is  made  that,  inherent 
in  Equation  12  for  upstream  differencing  of  left  to  right  flow, 

(<►  =  (4^  >  then  using  Equation  11  plus  additional  rear- 

rangement  gives 
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When  0=0,  Equation  26  reduces  to  Equation  15,  as  expected.  This 
form  is  convenient  for  the  linear  stability  analysis  and  propagation 
tests  using  double-sweep,  tridiagonal  matrix  solution  procedures. 
Linear  stability  analysis  of  QUICKIST 

32.  Use  of  the  von  Neumann  method  resulted  in  the  following  ex¬ 
pression  for  the  amplification  factor 


A(a)  = 


D  exp  (Fla)  +  E  +  F  exp  (-Fla)  +  H  exp(-n2a) 

o  o 

A  exp(na)  +  B  +  Cl  exp  (-Fla) 


(27) 


A  simple  computer  program  written  in  complex  arithmetic  was  employed  to 
find  [a|  and  Q  over  30  values  of  a  from  Equation  27  without  further 
simplification. 

33.  Some  results  when  0  =  0.5  are  shown  in  Figure  10  for 
C  =  0.5  and  Figure  11  when  C  =  1.0  .  In  both  instances,  the  explicit 


Figure  10.  Amplitude  portrait  for  QUICKEST  versus  QUICKIST  schemes 
when  C  =  0.5  ,  r  =  0  .  Both  have  no  phase  errors  for  all  of  or 
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Figure  11.  Amplitude  and  phase  portraits  for  QUICKEST  and 
QUICKIST  schemes  when  C  =  1,0  ,  F  =  0 

scheme  is  more  accurate  for  pure  advection.  The  extent  of  numerical 
damping  when  C  =  1.0  for  the  QUICKIST  algorithm  can  be  considered  ex 
cessive.  However,  as  shown  in  Figure  4  for  C=0.5,  r=l.l, 
QUICKIST  (0  =  0.5)  appears  more  accurate  for  all  a  .  Figure  5  also 
reveals  that  when  C  is  again  0.5  but  F  now  1.2,  QUICKIST  is  stable 
whereas  the  explicit  scheme  is  unstable.  The  implicit  algorithm  has 
far  superior  phase  response  for  these  latter  examples.  Finally,  Fig* 
ure  12  reveals  that  the  explicit  scheme  is  more  accurate  over  all  a 
for  the  pure  diffusion  case  when  F  =  1/6  . 

34.  The  (F,C)  plane  was  mapped  for  the  implicit  scheme  to 
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Figure  12.  Amplitude  response  for  pure  diffusion  (C  =  0)  when 
r  =  1/6  for  both  QUICKEST  and  QUICKIST  schemes 

determine  its  stability  range*  for  comparison  with  Leonard's  results 
in  Figure  3.  The  results  were  somewhat  disappointing.  It  was  antici¬ 
pated  that  the  implicit  QUICKIST  algorithm  would  have  a  much  larger 
stable  region  than  shown  in  Figure  13.  However,  the  unstable  corner 
below  C  =  0.5  and  beyond  P  =0.5  was  eliminated.  The  stable  range 
of  QUICKEST  method  is  also  shown  in  Figure  13  for  comparison. 

35.  The  amplitude  accuracy  of  QUICKIST  was  characterized  in  the 
stable  range  out  to  P  =  1.3  using  the  RMS  difference  norm,  as  before. 
Figure  14  presents  the  "contour"  map  of  RMS  difference  values  and  gives 
a  weighted-average  value  in  the  stable  area  (arbitrarily,  P  <  1.3)  of 
0.0946.  As  expected,  QUICKIST  with  0  =  0.5  is  less  accurate  than  the 
explicit  Leonard  algorithm. 


*  |A|  >  1.001  for  any  0  <  a  <  71  was  the  criteria  employed. 
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Figure  13.  Stable  range  of  QUICKIST  method  in  (r,C)  plane 
as  determined  in  this  study.  Also  shown  is  stable  range 
of  QUICKEST  scheme  and  the  QUICK8ST  algorithm  described 
later  in  this  report 
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Figure  14.  Iso-amplitude  lines  of  equal  RMS  difference  values  in 
(r,C)  plane  for  QUICKIST.  Stability  range  also  shown 


Propagation  tests  of  QUICKIST 

36.  The  QUICKIST  inaccuracies  are  confirmed  by  results  of  the 
triangle  propagation  tests.  Excessive  numerical  damping  results  to  the 
extent  demonstrated  in  Figure  15  for  C  =  0.5  (Fig.  15a)  and  C  =  1.0 
(Fig.  15b),  When  compared  with  Figure  7  and  the  fact  that  C  =  1.0 
gives  an  exact  result  for  the  explicit  scheme,  it  must  be  concluded  that 
the  implicit  QUICKIST  algorithm  is  too  inaccurate  to  warrant  further 
study. 


(b)  When  C=1.0  and  after  100  time  steps. 

Figure  15.  Propagation  tests  for  pure  advection  using 
periodic  boundary  conditions  and  QUICKIST  algorithm. 
Large  numerical  diffusion  present 
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Fully  Centered,  Implicit  Algorithm  (QUICKeST) 


37.  The  inaccuracy  of  the  QUICKIST  algorithm  can  be  traced  to  the 
off-centering  of  the  derivative  at  (i,n).  The  sole  advantage  of  this 
method  is  that  when  0=0,  the  explicit  QUICKEST  algorithm  is  regained. 

38.  A  fully  centered,  implicit  method  would  utilize  grid  point 
(i,n  +  1/2)  in  the  Taylor  series  expansion  so  that  even  the  0^  deriva¬ 
tive  is  effectively  weighted  by  0  .  When  0=0,  the  explicit  QUICK¬ 
EST  scheme  would  no  longer  be  recovered.  However,  far  greater  numerical 
accuracy  is  anticipated  since  all  terms  and  truncation  errors  are  cen¬ 
tered  at  the  same  time  level  by  0  .  Consequently,  this  fully  centered 
implicit  Leonard  scheme  is  given  the  acronym  QUICK0ST. 

Derivation 

39.  Consider  the  schematic  grid  sketched  in  Figure  16.  A  two- 
dimensional  Taylor  expansion  of  the  form 


f(a  +  h,  b  +  k)  =  f(a,b) 


'  '  x=a 

t=b 


+ 


) 


f (x,t) 


+  . . . 


x=a 

t=b 


is  utilized  where 


a  =  i 

b  =  n  +  1/2 
h  =  Ax 
k  =  At/2 


(28) 


(29) 


for  the  nomenclature  in  Figure  16.  The  initial  expansion^*  about 
(i,n  +  1/2)  yields  by  analogy  with  Equation  23 


*  A  coefficient  was  discovered  to  be  missing  in  this  expansion  at  a  later 
date.  The  complete  results  for  the  correct  expansion  are  given  in 
Appendix  A.  Surprisingly,  the  correct  expansion  (QUICK0ST  2  exhibited 
little  if  any  improvement  over  the  original  explicit  scheme  QUICKEST. 
Therefore,  it  was  dropped  in  favor  of  the  "chance*'  algorithm  QU1CK0ST 
discussed  in  detail  herein.  It  was  not  possible  to  determine  the  cause 
for  this  unexpected  result  within  the  allotted  scope  of  the  study. 
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Figure  16.  Schematic  of  operator  for  QUICK6ST 
scheme  determination 
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The  truncation  error  terms  0,  0,  0,  and  @  are  now  centered 

at  (i,  n  +  1/2)  giving  different  coefficients  and  incorporating  mixed 
time/space  derivatives.  Neglecting  all  fourth  and  higher  order  deriva¬ 
tives;  using  Equations  6,  7,  and  8  to  convert  all  time  and  mixed-time 
derivatives  to  pure  space  derivatives;  and  combining  all  like  deriva¬ 
tives  for  truncation  errors  yields 


.n+1  n 
(J).  -  4>. 

At 


+  u(l  -  6) 


.  „e  i"  -  <-! 


2Ax 


2Ax 


K(1  -  6)Ji3 - ! - LJ.  t  k6  - llll 


Ax 


n+1/2 


0  +  0 

t  (I  -  26)] 

@  +  0  +  0  +  0+  0 

^  fuAx^  5Ku(l  -  2e)AtlA 

^  te  12 - ^J('*’xxx). 


Ax" 


n+1/2 


(31) 


In  Equation  31,  truncation  error  from  terms(^and(^contribute  to  ^ 
while  all  truncation  error  terms  contribute  to  the 


^xxx 


XX 

coefficient. 
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40.  It  is  now  assumed  that: 

a.  The  truncation  errors  in  Equation  31  are  evaluated  at 
time  level  n  . 

b.  The  term  is  centered  upstream  (grid  point 

(i  -  1/2)  for  left  to  right  flow). 

Assumption  a  is  required  for  I'O  retain  the  tridiagonal  structure 

of  the  algorithm  as  in  QUICKIST.  Also,  it  is  felt  that  for  consistency, 
should  also  be  taken  at  level  n  .  Assumption  b  follows  from 
Part  II  to  produce  the  quadratic,  higher  order,  accurate,  upstream  dif¬ 
ference  operator  for  advection.  Using  these  two  assumptions  and  making 
the  usual  rearrangements  gives 


+  E<t>i  + 


where 

A  =  f  -  re 
B  =  1  +  are 
Cl  =  -  (^  +  re) 

2 

D  =  [-  §  (1  -  e)  +  r(i  -  0)  +  ^  (1  -  ae) 

*5-12  ‘='■(1  -  26)] 

E  =  [i  -  2r(i  -  e)  -  S-  (1  -  26)  - 1 1 1|  crd  -  26)] 
F  =  [I  (1  -  0)  r(i  -  6)  +  ^  (1  -  ae) 

+  I  -  ^  cr(i  -  ae)] 

"  =  -  [5  ■  T2  '='■(>  -  29)] 


(3a) 


(a) 

(b) 

(c) 


(d) 

(e) 


(f) 

(g) 


The  coefficients  D  ,  E  ,  F  ,  and  H  in  Equation  3a  are  different 
from  the  first  implicit  scheme  given  by  Equation  a6.  As  an  example, 
consider  the  coefficients  when  C=1.0,  r=0,  0=0. 5 


QUICK0ST 


1  .n-H  aH+I  1  ^  A®  +  6  A*'  +  2_  A®  _  ^  A® 

4  '•'i+1  ^  ♦i  '4  ♦i-i  ~  ■  la  ^i+i  ^  la  ^  ^  Ta  Vi  12  h-2 
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QUICKIST 


1  ^n+l  .  aD+I  1  ,n+l  1  iD  .  3  . n 
4  Vl  ^  '►i  '4  ♦i-l  =  4  4  ^i-l 


QUICKEST  (6  =  0) 
♦“*'  -  ♦“  , 


(b)  (33) 


(c) 


It  appears  that  QU1CK6ST  retains  more  of  the  intended  upstream  interpo¬ 
lation  flavor  of  the  Leonard  QUICK  method. 

Linear  stability  analysis  of  QUICK0ST 

41.  Equation  27  is  again  employed  but  with  coefficients  calcu¬ 
lated  from  Equation  32.  Some  results  (0  =  0-5)  are  shown  in  Figure  17. 
This  implicit  scheme  is  exact  in  amplitude  response  for  all  a  when 

C  =  1.0  or  0.5  for  pure  advection-  Some  phase  errors  exist  for 
higher  wave  numbers  as  demonstrated  in  Figure  17(b).  Figure  17  when 
compared  with  Figures  10  and  11  reveals  that  QUICK6ST  is  vastly  superior 
to  QUICKIST  and  similar  to  the  explicit  scheme  for  these  Courant  num¬ 
bers.  Both  implicit  schemes  reduce  to  the  Crank-Nicholson  scheme  for 
pure  diffusion  (C  =  0,  0  =  0.5). 

42.  Mapping  the  stable  area  of  QUICK0ST  in  the  (r,C)  plane 
revealed  a  much  larger  and  uniformly  stable  region  as  shown  previously 
in  Figure  13.  QUICK0ST  is  stable  for  all  P  below  C  =  1,5  except 
in  the  small  area  above  C  =  1.0  near  the  ordinate. 

43.  The  amplitude  accuracy  of  QUICK0ST  is  summarized  in  Figure  18 
out  to  r  =  1.3  again  using  the  RMS  difference  as  norm.  The  weighted- 
average  value  in  the  enlarged  stable  area  (P  <1.3)  is 


QUICK0ST 

0.0890 

and,  for  comparison, 

QUICKIST 

0.0946 

QUICKEST 

0.0568 

As  expected,  the  QUICK0ST  algorithm  is  more  accurate  in  amplitude 
response  over  all  wave  numbers  than  QUICKIST.  It  is  generally  more 
inaccurate  than  the  explicit  scheme,  but  is  far  more  useful  due  to 
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a 


(b)  Some  phase  errors  present. 

Figure  17.  Amplitude  and  phase  portraits  for  QUICK6ST  scheme  and 
pure  advection  (F  =  0)  case  when  C  =  0.5  and  1.0 


r 

Figure  18.  Iso-amplitude  plot  for  QUICK6ST  with  stable  range 

also  shown 
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the  extended  stable  region  above  C  =  1.0  for  practical  values  of 
r  >  0.5  .  The  phase  accuracy  for  the  three  schemes  generally  followed 
these  same  trends. 

Propagation  tests  of  QUICK0ST 

44.  The  excellent  numerical  accuracy  of  the  QUICK0ST  is  con¬ 
firmed  by  the  triangle  propagation  test  results  depicted  in  Figure  19. 

The  discrepancy  at  the  peak  is  comparable  to  that  produced  by  the  ex¬ 
plicit  scheme  (Figure  7)  and  the  waviness  of  the  upstream  face  is  due 
to  phase  errors  (recall  Figure  17),  The  closeness  of  numerical  results 
between  the  explicit  and  QUICK6ST  algorithms  is  also  demonstrated  in 
Figures  20  and  21  where  some  physical  diffusion  (F  =  0.1)  is  now  present. 


(b)  When  OKO  and  after  100  time  steps. 

Figure  19.  Propagation  test  results  for  pure  advection 
case  using  periodic  boundary  conditions  and  QUICK6ST 

algorithm 
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(a)  When  C=0.5,  r=0. 1,  P^=5  and  after  200  time  steps. 


(b)  When  C=1.0,  r=0. 1,  and  after  100  time  steps. 


Figure  20.  Propagation  test  results  for  example  with  some  physical 
diffusion  present  (T  =  0.1)  comparing  explicit  (dots)  and  implicit 
(triangles)  QUICK6ST  algorithms 


In  these  examples,  only  initial  conditions  are  shown. 

45.  Other  propagation  tests  gave  similar  results. 
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Initial  Conditions 
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Figure  21.  Propagation  test  results  for  narrow  triangle  (JW  =  2)  and 
physical  diffusion  (F  =  0.1)  comparing  QUICKEST  (dots)  and  QUICK6ST 
(triangles)  algorithm  when  C=0.5,  r=0.1,  after  200 

time  steps 


PART  IV:  CONCLUSIONS  AND  RECOMMENDATIONS 


46.  The  purpose  of  this  project  was  to  develop  and  study  numeri¬ 
cal  properties  of  an  implicit  algorithm  consistent  with  the  QUICKEST 
scheme  of  Leonard  (1979)  for  the  unsteady,  one-dimensional  transport- 
diffusion  equation.  Based  upon  the  results  discussed  in  this  report, 
the  following  conclusions  are  drawn  and  recommendations  offered. 

Conclusions 


47.  The  explicit  algorithm  labeled  QUICKEST  is  identical  to  a 
forward-time,  centered-space  scheme  with  all  truncation  errors  ex¬ 
pressed  as  spatial  derivatives  and  removed  through  0  and  (|) 

XX  XXX 

terms.  The  term  is  centered  upstream  and  results  in  a  higher 

order,  accurate,  upstream  difference  method. 

48.  A  fully  centered,  implicit  scheme  called  QUICK0ST,  as  devel¬ 
oped  herein,  is  of  comparable  accuracy  to  the  explicit  operator  yet  ex¬ 
tends  the  stable  range  in  the  (r,C)  plane  to  essentially  below  the 
line,  C  =  1.5  for  all  P  . 

49.  QUICK8ST  is  a  stable,  accurate,  and  robust  algorithm  of  con¬ 
siderable  engineering  usefulness  for  transport-diffusion  computations. 

50.  Because  development  of  the  implicit  QUICK0ST  algorithm  has 
proven  to  be  highly  successful,  further  research  on  the  QUICK6ST  2 
scheme  is  not  warranted  except  to  learn  the  reasons  for  its  poor 
performance . 


Recommendations 

51.  The  QUICK6ST  scheme  for  one-dimensional  flows  should  be  ex¬ 
tended  to  variable  grid  spacing.  Only  nomenclature  and  programming  mod¬ 
ifications  are  required. 

52.  A  two-dimensional  QUICK0ST  algorithm  should  be  developed. 

Use  of  an  alternating-direction  implicit  (ADI)  structure  should  be  con¬ 
sidered  along  with  the  **double-sweep’*  solution  algorithm. 
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53.  An  experimental  test  program  must  be  devised  to  fully  verify 
the  two-dimensional  QUICK6ST  against  existing  and  new  data  from  both 
the  field  and  the  laboratory. 

54.  Systems-oriented  computer  codes  and  extensive  software  should 
be  written  for  ease  in  application  by  computational  hydraulic  engineers. 
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APPENDIX  A:  DERIVATION  OF  QUICKOST  2 


1.  As  mentioned  in  the  footnote  on  page  37,  the  original  expansion 
to  develop  the  QUICK0ST  algorithm  discussed  in  this  report  was  found  to 
be  incorrect  (Ross  Hall,  19  Oct,  1982,  personal  connminication) .  This 
Appendix  gives  the  explanation,  all  corrected  equations,  and  a  linear 
stability  analysis  of  the  results.  This  corrected  version  is  dubbed 
QUICK0ST  2  and  is  shown  to  be  a  far  less  robust  algorithm  than  the 
*'chance*'  algorithm  QUICK0ST.  For  this  reason,  the  original  report  and 
all  conclusions  were  felt  to  still  be  valid  and  the  corrected  algorithm 
QUICK0ST  2  relegated  to  this  Appendix, 

2.  Consider  the  schematic  operator  shown  in  Figure  16,  The 
general  two-dimensional  Taylor  series  expansion  given  in  Equation  28  is 
correct  but  the  indices  should  be 

a  =  i 

b  =  n+0 

h  =  ±AX 

k  =  +(l-0)At,-0At  (Al) 

This  will  give,  for  example,  a  space  difference  from  n+0  to  n+1  (over 
+(l-0)At)  and  from  i  to  i+1  (over  +Ax) 

n+1  ,n+0  ,  r*  ^n+0  m  \n+0n 

^i+1  '  lAx{4>^)^  +  {l-e)At(4>^)^  J 


Al 


The  coefficients  underlined  in  Equation  A2,  namely  the  ^  and  were 
inadvertantly  left  out  of  this  and  all  other  expansions  in  the  original 
derivation  that  produced  the  **chance”  algorithm  Q13ICK6ST. 

3*  The  correct  expansion  gives  the  following  equations. 
Corrected  Equation  30  (pp  3B-39) 


(1) 


(3) 
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Corrected  Equation  31  (p  39) 
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=  K(1-0)-1^3 - } - L_L  +  Ke-lll 
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(l-20)uKAt](<(>^^^) 
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where  Term  (8)  +  Term  (10)  _  0 

Corrected  Coefficients  in  Equation  32  (p  40) 

A  =  -  r0-  R0 

B  =  1  +  2r0  +  2R0 
Cl  =  -(^  +  r0  +  R0) 

D  =  ---^-2-—  +  r(1-0)  +  R(l-0)  +  s 
E  =  1  -  2r(l-0)  -  2R(l-0)  -  3S 
F  =  +^^^^  +  r(l-0)  +  R(l-0)  +  35 

H  =  -S  (A5) 

where : 

R  =  ^  (1-28) 

s  =  C’  -  (l-2e)Cr]  (a5) 


Corrected  Equation  33a  (p  40)  when  C=l.O,r=0,  0=0.5,  QUICK0ST  2 


1  n+1 

4  ‘*’1+1 


C’ 


1  ,n+l 


4  ‘*’1-1 


+  t) 


n  1 
i-1  ■  ?  *1-2 


(A6) 


4.  A  linear  stability  analysis  of  QUICK0ST  2  using  Equation  27 
was  employed  but  with  coefficients  calculated  from  Equation  A5  above. 
Some  results  are  summarized  below  for  a  perfectly  centered  scheme  with 


0=0.5, 

c 

r 

Remarks 

1.0 

0 

Exact  for  all 

0.75 

0 

Stable 

0.5 

0 

Stable 

1.0 

1.0 

Stable 

A3 

c 

T_ 

Remarks 

1.01 

1.0 

Unstable 

1.5 

0 

Unstable 

2.0 

0.5 

Unstable 

5.  The  scheme  QUICKGST  2  becomes  unstable  for  all  T  when  the 
Courant  number  exceeds  unity  which  is  an  unexpected  result  for  an 
implicit  scheme.  Using  other  values  of  0  (0.4,  0.51,  0.75,  1.0) 
produced  erratic  results  with  a  stable  region  always  far  smaller  than 
that  given  by  the  QUICK6ST  scheme.  We  conclude  that  this  algorithm  is 
far  less  robust  and  inferior  to  the  QUICK0ST  algorithm  discussed  in 
this  report. 
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